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A FRACTIONAL SPACE-TIME OPTIMAL CONTROL PROBLEM: 
ANALYSIS AND DISCRETIZATION* 

HARBIR ANTILt, ENRIQUE OTArOLA*, AND ABNER J. SALGADO§ 

Abstract. We study a linear-quadratic optimal control problem involving a parabolic equation 
with fractional diffusion and Caputo fractional time derivative of orders s G (0,1) and 7 G (0,1], 
respectively. The spatial fractional diffusion is realized as the Dirichlet-to-Neumann map for a 
nonuniformly elliptic operator. Thus, we consider an equivalent formulation with a quasi-stationary 
elliptic problem with a dynamic boundary condition as state equation. The rapid decay of the 
solution to this problem suggests a truncation that is suitable for numerical approximation. We 
consider a fully-discrete scheme: piecewise constant functions for the control and, for the state, first- 
degree tensor product finite elements in space and a finite difference discretization in time. We show 
convergence of this scheme and, for s G (0,1) and 7 = 1 , we derive a priori error estimates. 

Key words, linear-quadratic optimal control problem, fractional derivatives and integrals, 
fractional diffusion, weighted Sobolev spaces, finite elements, stability, fully-discrete methods. 
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1. Introduction. We are interested in the design and analysis of efficient so¬ 
lution techniques for a linear-quadratic optimal control problem involving an initial 
boundary value problem for a space-time fractional parabolic equation. Let 12 be 
an open and bounded domain in M" (n > 1), with boundary (9f2. Given s € (0,1), 
7 € (0,1], and a desired state : f2 x (0, T) —>• M, we dehne 

= \J^ (llu - Ud||i2(o) +f^Mhin)) dL 

where /i > 0 is the so-called regularization parameter. Let f : 12 x (0, T) —>■ M and 
uq : n —t M be fixed functions. We will call them the right hand side and initial datum, 
respectively. We shall be concerned with the following optimal control problem: Find 

minJ(u,z), (1-1) 

subject to the space-time fractional state equation 

97u -I- r®u = f -f z in f2 X (0, T), u(0) = Uq in 12, (1.2) 

and the control constraints 

a{x',t) < z{x',t) < b{x',t) a.e. {x',t) G Q := ft x {0,T). (1.3) 

The functions a and b both belong to L‘^{Q) and satisfy the property a (a;', t) < b{x', t) 
for almost every [x',f) G Q. The operator with s € (0,1), is the fractional power 
of the second order elliptic operator 

Cw = —diva;/(AVa;''u;) -|- cw in 12, w = 0 on (912, (1.4) 
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where 0 < c G and A G GL(n,M)) is symmetric and positive definite. 

The fractional derivative in time for 7 G (0,1) is understood as the left-sided 
Caputo fractional derivative of order 7 with respect to t, which is formally defined by 


dfu{x\t) = 


1 


du{x', r) 


r(l-7) Jo (t-rp 


dr 


dr, 


(1.5) 


where T is the Gamma function. For 7 = 1 , we consider the usual derivative dt- 


For convenience, we will refer to the optimal control problem (1.1)-(1.3) as the 


space-time fractional optimal control problem] see section for its precise description 
and analysis. One of the main difficulties in the study of the state equation (1.2) is 


the nonlocality of the fractional time derivative and the fractional space operator (see 
01711111131133]). A possible approach to overcome the nonlocality in space is given by 
the result of Caffarelli and Silvestre in M” 0 and its extensions to bounded domains 
0 133] : Fractional powers of the spatial operator C can be realized as an operator 
that maps a Dirichlet boundary condition to a Neumann condition via an extension 
problem on the semi-infinite cylinder C = 11 x (0,(X)). Therefore, we shall use the 


Caffarelli-Silvestre extension to rewrite the fractional space-time state equation (1.2) 
as a quasi-stationary elliptic problem with a dynamic boundary condition: 

J-div(y“AV'i') = 0 inC x ( 0 ,r), = 0 on drC x ( 0 ,r), 

ydf^ -\- J^df^ = f -|- z on (12 X {0}) x (0, T), = Uq on 12 x {0}, t = 0, ^ ^ 

where drC = 912 x [0,oo) is the lateral boundary of C, a = 1 — 2s G (—1,1), ds = 
2 “r(l — s)/r(s) and the conormal exterior derivative of at 12 x { 0 } is 


df^ = - lim y°‘%, 


(1.7) 


where the limit must be understood in the distributional sense 00133]. Finally, 
A{x',y) = diag{A(a;'), 1} G GL(n-|-1, K)). We will call y the extended variable 

and the dimension n -I- 1 in the extended dimension of problem ( |1.6[ ). As noted 

in 00133], and the Dirichlet-to-Neumann operator of (|1.6l) are related by 


dsC‘'u = 9“^ in (12 x {0}) x (0, T). 


We briefly elaborate on these ideas in (2.5 A rigorous analysis is provided in mm- 
The study of solution techniques for elliptic and parabolic problems involving 
fractional derivatives is a relatively new but rapidly growing area of research. We 
refer the reader to 11711131130] for an overview of the state of the art. Numerical 
strategies for solving a discrete optimal control problem with PDF constraints have 
been widely studied in the literature; see [n na m US] for an extensive list of 
references. Mainly, these references are concerned with control problems governed by 
elliptic and parabolic PDFs, both linear and semilinear. The common feature here is 


that, in contrast to (1.1|-(1.3), the state equation is local. 


The numerical analysis of optimal control problems involving evolution equa¬ 
tions with fractional diffusion and fractional time derivative is still at its infancy. 
To the best of our knowledge, the first work that provides a comprehensive treat¬ 
ment of an optimal control problem involving fractional elliptic operators in space is 
0. Concerning fractional derivatives in time, the first work that attempts to study 
an optimization problem constrained by a fractional order ODF is 0 where, through 
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completely formal calculations, the author derives optimality conditions and a numer¬ 
ical scheme. However, no justification is provided for either the optimality conditions 
nor the numerical scheme. Later, similar optimization problems have been discretized 
via a finite element method [5], a modified Griinwald-Letnikov approach IS [ID] and 
a rational approximation approach |34j . However, fundamental mathematical results 
such as stability and convergence of the proposed numerical schemes are missing in 
these works. Recently, convergence of spectral based techniques has been explored 
in [5S US for an optimization problem restricted to fractional order ODEs. Opti¬ 
mal control problems for one dimensional evolution equations with only fractional 
time derivatives have been recently studied in ESEH]. In these references, the au¬ 
thors derive rigorously first order necessary optimality conditions, propose numerical 
schemes based on spectral methods and obtain a priori error estimates. These error 
estimates, however, are derived using regularity assumptions that are verified only in 
very restricted cases muzi. 

We provide a comprehensive treatment of a linear-quadratic optimal control prob¬ 
lem involving evolution equations with fractional diffusion and fractional time deriva¬ 
tive: s G (0,1) and 7 S (0,1]. To the best of our knowledge this is the first work 
addressing such a problem from a mathematical point of view. We rigorously derive 
optimality conditions, present a numerical scheme and prove its convergence. In ad¬ 
dition, for s G (0,1) and 7 = 1 , we derive a priori error estimates. We overcome 
the nonlocality of by using the results of Caffarelli and Silvestre [ 6 ]. We realize 
the state equation by ( | 1 . 6 [ ) so that, our problem can be equivalently written as: 
Minimize J subject to the extended state equation (1.6) and the control constraints 


(1.3). 


Inspired by 0123 US, we propose a simple strategy to find the solution to the 
space-time fractional optimal control problem (1.1)-(1.3): given f and u^, we realize 


(1.2) by (1.6) and apply standard techniques to solve this problem. We thus obtain an 
optimal control z : H x (0, T) —K and an optimal state ^ : C x (0, T) —>■ K. Letting 
u : H X (0,T) 9 (a:',t) ‘^{x',0,t) e K we obtain (u,z) that solves (1.1)-(1.3). 

The outline of this paper is as follows. In section 0 we introduce notation, recall 
elements from fractional calculus, define fractional powers of elliptic operators via 
spectral theory and show the equivalence with the Caffarelli-Silvestre extension. This 
allows us to study (1.6) and provide some energy estimates. On the basis of this, 


in section 0 we study the space-time fractional optimal control problem. We derive 
existence and uniqueness results together with first order sufficient and necessary 
optimality conditions. In Q we begin the numerical analysis of our problem. We 
introduce a truncation of the state equation and derive approximation properties of 
its solution. In section]^ we recall the fully discrete scheme of [27] that approximates 
the solution to the state equ ation (1.2). For s G (0,1) and 7 = 1, we derive a novel 
L^(Q)-error estimate in ^5.4] Subsection 6.1 is devoted to the design of a numerical 


scheme to approximate the control problem IPll-dOt, and in | 6 . 2 [ we derive a priori 
error estimates for s G (0,1) and 7 = 1 . The convergence of the scheme is analyzed 
in ]|6.3|for s G (0,1) and 7 G (0,1]. Finally, sectionpresents numerical experiments 


that illustrate the theory developed in 

2. Notation and preliminaries. Let us set notation and recall some facts that 
will be useful later. 


2.1. Notation. Throughout this work fl is an open, bounded and connected 
subset of K", n > 1, with polyhedral boundary 917. If T > 0 is a fixed time, we set 
Q = 17 X (0,T). We will follow the notation of [2Z1 [221 and define the semi-infinite 
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cylinder with base and its lateral boundary, respectively, by C = x (0, oo) and 
SlC = dVl X [0, oo). For 9^ > 0, we define the truncated cylinder Cy = Vlx (0, fK) and 
dyCy accordingly. Since we will be dealing with objects defined on M" and it 

will be convenient to distinguish the extended n+ 1-dimension. If a; S we write 

X = (a;', y), with x' S K" and y S M. 

If X and y are normed spaces, X ^ y means that X is continuously embedded 
in y. We denote by X' and II-lU the dual and norm of X ^ respectively. The relation 
a b indicates that a < Cb, with a nonessential constant C that might change at 
each occurrence. 

If £> C is open, fV > 1, and : U x (0, T) —)■ M, we will regard (j) as a 
function of t with values in a Banach space X, i.e., : (0, T) 9 t >->• t) S X. 

For 1 < p < oo, LP{0,T; X) is the space of df-valued functions whose df-norm is in 
LP{0,T). This is a Banach space for the norm 

II0I|lp(o.T;A') = ( / ||(?i'(t)||^dt ] , l<p<oo, ||()'||Loo( 0 ,r;A') = esssup||0(t)|U. 

\Jo / te(o,r) 


2.2. Fractional derivatives and integrals. The left Caputo fractional deriva¬ 


tive is defined in (1.5). The right-sided Caputo fractional derivative is mm- 


dy-tgit) = - 


1 


_ r g'(o 


d^, 7 6(0,1). 


( 2 . 1 ) 


For g 6 L^(0,T) and a ^ 0, the cLird HlG-TTLOiTlTh L'IOU'VzIIg J'T'QiCi/IOTKxl ZTltGQT'Qils 
of order a are, respectively, [351 Definition 2.1, §2] 

[32l §2.2-2.3] provides a motivation for these definitions inspired by the Abel equation. 

Proposition 2.1 (continuity of fractional integrals). For a > 0 and 1 < p < oo, 
then If and are continuous from U‘{Q,T) into itself and 




ll^f-t5llLP(0,T) < _j_ 2^) II5||lp(0,T)- 


for all g 6 LP(0,T). These maps also are continuous from C([0,r]) into itself. 

Proof. For the proof of continuity in T^(0, T) see [311 Theorem 2.6, §2]. To obtain 
the continuity in C([0,r]) we use the continuity in L°°{0,T), together with the fact 
that if p 6 C'([0,T]) then its fractional integrals are continuous as well. This can be 
easily shown by recalling that g is also uniformly continuous. □ 

We also define the left and right Riemann-Liouville fractional derivatives of order 
7 6 (0,1), respectively, by |32l Definition 2.2, §2.3] 


D7g{t) 


1 d 
F(1 — 7 ) dt 


r gio 

Jo 


dC, D'fy_^g{f) 


-1 d 
F(1 — 7 ) dt 



g{0 

(e - tv 


df. 


A relation between the Caputo and Riemann-Liouville derivatives is given below. 

Lemma 2.2 (relation between fractional derivatives). Let 7 6 (0,1) and g 6 
Wf{0,T), then Dfg and Dy_^g exist almost everywhere on [0,T]. In addition, 
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Dtg,D'^-t9 S L^{0,T) /or 1 < r < i and 


D]g{t) = d^g{t) 


1 g(o) 

r(i- 7 ) ’ 


Dl_tgit) = + 


1 gjT) 

r(i- 7 ) {T-ty 


:• (2.3) 


Proof. See [32l Lemma 2.2, §2.3]. □ 

We now derive an integration by parts formula for Caputo derivatives that will 
be fundamental in our analysis. For 7 € (0,1) we define 

= {/ e C([0,T]) : 57 / G L'(0,r)}, R^ = {g€ C([0,r]) : d^_,g G L'(0,r)}. 

Lemma 2.3 (fractional integration by parts formula). /// G L.y and g G 
then the following fractional integration hy parts holds: 

r dyf{t)g{t) dt + /( 0 )( 4 i^ 5 )( 0 ) = r fit)df_,g{t) dt + giT){ll-y)iT). (2.4) 
^0 Jo 

Proof. If / and g are smooth, recall that [521 Corollary 2, §2.6]: 

f Dy{t)g{t)dt= f f{t)Df._^g{t)dt; 


(|2.4|) now follows from (|2.3|) and (|2.2|); the point values are well defined since / G 

Conclude by density. □ 


2.1 


C([0,T]) implies f,iy_^f G C([0,r]); see Proposition 

It is important to remark that there is another definition, not completely equiva¬ 
lent, of fractional derivatives: the so-called Griinwald-Letnikov derivative [TB]. Among 
all possible definitions of fractional derivatives, we adopt the left-sided Caputo frac¬ 
tional derivative as in problem (1.2): the Caputo approach leads to an initial con¬ 
dition of the form u = Uq which is physically meaningful. The Riemann-Liouville ap¬ 
proach leads to initial conditions containing the limit values of the Riemann-Liouville 
fractional derivatives at t = 0 , something that does not have a clear physical meaning. 

2.3. Fractional powers of second order elliptic operators. Spectral theory 
for the operator C yields the existence of {{Xk, <fk)}kGn C K"*" x Hq{Q) such that 

C(pk = Xk<Pk in ^2, tpk = 0 on 511, fc G N. (2-5) 

{‘/SfcjfeGN is an orthonormal basis of L^{fl). Fractional powers of £, are defined by 


OO p 

C‘^w:='^Xlwkipk, ywGC^iXl), sG(0, 1), Wk = 

k=i 


Wifk. 


( 2 . 6 ) 


By density, (2.6) can be extended to IHI®(11) = [L^(ll), ILq( n)]^. If we denote by 
IHI“®(11) the dual of IHI^(n), then : ]HI®(f2) —>■ IHI“^(n) is an isomorphism. 


2.4. Weighted Sobolev spaces. To study ( |1.6| ) we consider Sobolev spaces 
with the weight \y\°‘, a G (—1, 1). For D C we define 




with norm 


\W\\m(\y\a.^D) = 


D) = {w& y{\y\^,D) : |Vu;| G LWy^y)} 

(ii 




IVwf 


L^(\V\-‘.D) 


1/2 


(2.7) 
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Since a e (—1,1), |j/|“ belongs to the Muckenhoupt class see O |36]. 

Then, H^{\y\°‘,D) is Hilbert and C°°{D)nH^{\y\°‘, D) is dense in D) (cf. [Ml 

Proposition 2.1.2, Corollary 2.1.6], m and mi Theorem Ij). 

We also define the weighted Sobolev space 

= {wG H\y^,C) : zc = 0 on d^C) . (2.8) 

As |29[ (2.21)1 shows, the following weighted Poincare inequality holds: 

lk||L2(y._C) < l|Vw||L2(y._c), Vw e Hi(?/“,C), (2.9) 


thus ||Vzc||l 2 (j^c c) is equivalent to (2.7) in H\{y°',C). For w S H^{y°‘,C), ti^w 


denotes its trace onto H x {0}. We recall (UHl Prop. 2.5]) 

trni?i(y“,C) =H®(H), jj trn wjlHocn) < lkll^i(j^c_c)- 


( 2 . 10 ) 


2.5. The state equation. We follow m and define 

W:={wG L^{0,T;L'^{n)) n L 2 (o, T;H"(ll)) : d^w € L2(0, T;H-®(H))}, 

W:={wG L2(0,rAi(y“,C)) : t^w G L'^{Q,T;E.-^{n))}. 


The Caffarelli-Silvestre extension result for problem (1.2) reads 013 [Ml HZ] : Given 
f,z S L^(0,T;IHI“®(n)), the function u G W solves (1.2) if and only if its harmonic 
extension G V solves the following version of (1.6): Find 'W G W such that 
tro ^ (0) = Uo and for a.e. t G (0, T) 


(trn = (f + z,trn (/>) \/(j) G Hl{y°‘,C), 

where (•, •) is the duality pairing between IHI®(11) and IHI“®(H) and 

a{w,(j)) = \- I y^^ {A{x')Ww-Wcj) + c{x')w(j)) dx'dy. 
ds Jc 


( 2 . 11 ) 


( 2 . 12 ) 


The regularity of A and c implies that a is bounded and coercive in H]^{y°‘,C). In what 
follows, we shall use repeatedly that a{w,wy^^ is a norm equivalent to | • l^i(yo ^y 
Denote 


A2(v,g) := J^"||v|li.(^)(T) + ||g|li2(o.T;H-.(a)), (2-13) 

^^(v,g) := ||v|l|2(f2) + llg|li2(o,T;H-'’(n))> (2-14) 


where is the left fractional integral of order 1 — 7 defined in ( 2 . 2 ). 

Theorem 2.4 (existence and uniqueness of u and ^). Given s G (0, 1), 7 G (0, 1], 
f,z G L^(0,T;IHI“®(H)) and Uq G L'^{GL), problems ( |1.2[ ) and (2.11) have a unique 
solution. In addition, we have the following energy estimates for u, solution to 


It ^lhlli 2 (Q)(T) + ||u||| 2 (o,T;H'>(n)) ^ A^(uo,f + z), (2-15) 

ll^llL“( 0 .T;L 2 (a)) + l|u||i 2 (o,T;H»(n)) S^(uo,f+ z). (2-16) 


In addition, we have following energy estimates for ‘W solution to (2.11 ).• 


I troII 


L“(o,r;L2(n)) 


+ 


L^(0,T-Hl{y^,C))^^ (Uo,f + z), 


(2.17) 

(2.18) 
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where the hidden constants do not depend on u, ^ nor the problem data. 

Proof. The well-posedness of (1.2) and (1.6), together with (2.15) and (2.17) are 
presented in Theorem 2.6 and Corollary 3.8]. The estimates (2.16) and (2.18) 
follow from the arguments developed in Ham]. □ 

Remark 2.5 (7 = 1). Given g G LP{0,T), we have 1‘^g -g g in U'{Q^T) as cr 0 
[32l Theorem 2.6]. Take the limit as 7 f 1 in (2.15) and (2.17), to recover the well 
known energy estimates for parabolic equations with first order derivative in time. 

Remark 2.6 (continuity in time). An adaption of [31] Theorems 2.1-2.2] shows 
that, for every 7 e (0,1] and s G (0,1), the solution trn= u G C([0, T]; L^(n)). 
This is not only necessary to make sense of the initial condition, but also to derive 
optimality conditions, as we will see in section]^ 

We conclude with an elementary extension of Lemma |2.3| 

Lemma 2.7 (fractional integration by parts). Let 7 G (0,1]. // v,w G Wn 
C([0, r]; L^(r2)), then we have the following integration by parts formula: 


(97v(t), w(t)) - {^f_^w{t),v{t)) dt = (w(T), ^v)(T)) - (v(0), (4_»(0)), 


p-'y. 


where (•, •) = (•, •)l2(0) (•, •) is the duality pairing between IHI®(n) and El ®(n). 

Proof. When v and w are smooth we integrate (2.4). Conclude by density. □ 

3. The fractional control problem. In this section, we analyze the space- 
time fractional optimal control problem. We derive existence and uniqueness results 
together with first order necessary and sufficient optimality conditions. 

For J defined in (1.1) the fractional control problem reads: Find min J(u,z), 
subject to the state equation (1.2) and the control constraints (1.3). The set of 
admissible controls is defined by 

Zad := {w G LP{Q) : a{x',t) < w(x',<) < b{x',t), a.e. (x',<) G Q} , (3.1) 

which is a nonempty, bounded, closed and convex subset of L‘^{Q). To study this 
problem, following [331 §3], we introduce the control to state operator. 

Definition 3.1 (control to state operator). The map S : L^(0, T; IHI“®(n)) 9 z 1 —)■ 
u(z) G W, where u(z) solves (1.2) is called the fractional control to state operator. 


S is an affine and, by the estimates of Theorem |2.4[ continuous operator. More¬ 
over, since W ^ T^(Q) ^ T^(0,r;IHI“®(n)), we may also consider the operator S as 
acting from Lf‘{Q) into itself. For simplicity, we keep the notation S. We now define 
the optimal fractional state-control pair. 

Definition 3.2 (optimal fractional state-control pair). A state-control pair 
(u(z),z) G W X Zad is called optimal for the problem (1.1)-(1.3), ifb{z) = Sz and 

,/(u(z),z) < J(u(z),z), 


for all (u(z),z) G W x Zad such that u(z) = Sz. 

The existence and uniqueness of an optimal state-control pair is as follows. 
Theorem 3.3 (existence and uniqueness). The optimal control problem (1.1)- 
(1.3) has a unique solution (u(z),z) G W x Zad- 

Proof. Using the operator S, problem (1.1)-(1.3) reduces to: Minimize 

f{z):=hsz-Ud\\UQ, + ^\\z\\l,,QP ( 3 - 2 ) 


over Zad. Since p, > 0 the strict convexity of / is immediate. S is continuous, so / is 
weakly lower semicontinuous. Zad is weakly sequentially compact. The direct method 
of the calculus of variations (3] Theorem 1.15] allows us to conclude. □ 
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3.1. Formal Lagrangian formulation. We now formally derive first-order nec¬ 
essary and sufficient optimality conditions for the control problem dllll-dOl). We 
proceed via the Lagrangian approach described in §3.1]. We must emphasize 
that, although these computations are merely formal, they are quite insightful as 
they allow us to determine what is the correct form of the optimality conditions with 
a simple and intuitive procedure. 

Let p denote the adjoint variable, the Lagrangian L : W x Zad x W —>■ M is 

L(u, z, p) I-)-J(u, z) — f -b— f — z) p. 

JQ 

We expect the following necessary and sufficient optimality conditions [35l §3.1]: 


[ LipL(u,z,p)/i = 0 V/ieW, (3.3) 

JQ 

f ZljL(u, z, p)/i = 0 V/i S W, with/i(0) = 0, (3.4) 

JQ 

[ HzL(u,z, p)(z - z) > 0 VzeZad. (3.5) 

JQ 


We start with a formal computation which uses the integration by parts formula (2.4): 

[ {d^u + C^u-f -z)p= [ (9j_jp-l-r^p) u - f (f-l-z)p 
J Q J Q -J Q 

- [ u(o)(4i^p)(o)+ [ p(r)(//-"u)(r). 

Jn Jn 


Based on the previous computation, we rewrite expression (3.4) as follows: 


/ L>uL(u, z, p)h = (u - Ud, h)L^Q) - (^r-iP + h)L^{Q) 

JQ 


-(p(T),(/i-^M(T))^.(n)=0, 


(3.6) 


for all /i G W such that h{0) = 0. 

Let (j) G T), ijj G C^{n) be arbitrary and define h = where ip solves the 

Abel equation = (^(t). This is possible because of the unique solvability of 

the Abel equation given in [351 Theorem 2.1, §2.2]. Notice that with this definition 


[l] '^h){T) = 0. Using this particular choice of h in (3.6) yields 


^ ‘fit) {d^-tP + ^"P - (u - Ud), dt = 0. 


Owing to the results of [351 Theorem 13.2, Theorem 13.5], the range of the fractional 
integral ll~'^ contains all smooth functions. In other words the relation above must 
hold for all smooth and compactly supported which implies 


9J_(P -f £®p = u - Urf. 


(3.7) 


It remains to obtain a terminal condition for p. To do so, we notice that we have 

(p(T),(/^^/.)(T))^,(^)=0. 
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If we were allowed to set h constant in time this would yield p(T) = 0. However, since 
h{Q) = 0, the only admissible and constant in time function is h = Q. To circumvent 
this we set h = fe(t)x with y G (7^(0) arbitrary and i^{t) G C'°d([0,T]) given by 

4(i) = 0<t<eT, 4(t) = l, eT<t<T. 


This particular choice of h yields 

(p(T),x)^.(^) (//-^4)(T)=0. 

To conclude, it remains to notice that lime_,.o(/j^~^4)(7") = > 0- Collect¬ 

ing the derived equations, our formal argument yields the following strong system for 
the adjoint variable p. 

Definition 3.4 (fractional adjoint state). The solution p = p(z) G W o/ 


^T-iP + '^^P = u — in Q, p(T) = 0 in ft, 


(3.8) 


for z G L^(0,T;IHI ®(H)), is called the fractional adjoint state associated to u — u{z). 
Remark 3.5 (7 —>■ 1). For g G 1T^^(0,T), d'^_^g —)■ —dtg as 7 f 1 and we recover 

-dtP + £®p = u - Ud in Q, p(r) = 0 in H, 

a standard backwards parabolic problem with terminal condition. 

Well-posedness of (3.8) follows from a change of variables. If i/; : [0, T] —>• K define 
ffit) := ijj(T — t) and notice that i/’'(t) = —ij'iT — t). Therefore, with c-y = r(l — 7), 


rT-t 


Cyd2^{T-f) = 


[{T-t)-0^ 


de = - 


{g-ty 


dg = Cjd^_^'ip{t). 


As a consequence, the backwards in time problem ( |3.8[ ) with a right Caputo fractional 
derivative can be equivalently written as a forward in time problem with a left Caputo 


fractional derivative as (1.2|. The well-posedness of (3.8) then follows from (2.5 


We conclude this formal analysis with the following variational inequality: 

(/rz-f p,z-z)^ 2 (Q) >0, VzGZad, (3.9) 


which follows from ( |3.5| ). 

Remark 3.6 (Lagrangian approach). Although formal, this approach is system¬ 
atic and useful to derive optimality conditions of a control problem, specially in our 
case, where the state equation (1.2) involves fractional derivatives in time and space. 

3.2. Optimality conditions. We begin with a classical result. 

Lemma 3.7 (variational inequality). Let f be defined by (3.2). The function 
z G Zad minimizes the functional f if and only if 


(/'(z),z-z)l2(q) > 0, 


(3.10) 


for every z G Zad- 

Proof. See [35l Lemma 2.21]. □ 

To derive first-order optimality conditions, we need the following result. 

Lemma 3.8 (auxiliary result I). Let z denote the optimal control given by Theo¬ 
rem 


3.3 and u = Sz. Then, for every z G Zad, we have 


(U - Ud, U - u)i2(Q) = (p, Z - z)i2(Q) 


(3.11) 
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where u = Sz S W and p = p(z) S W solve problems (1.2) and (3.8), respectively. 
Proof. Define </> := u — u G W and notice that (j){0) = 0 in D. Moreover 


{d'f<i) + C(j3,w) = {z-z,w)Li(a) ywGW{n), a.e. (0,T). 
Since p G W setting w = p{t) in (3.12[) and integrating over time yields 


(3.12) 


(p,z-z)i2(Q) = /" (57())+ £*<(), p). 
^0 


Lemma 2.7 and the fact that the operator is self adjoint allow us to write 

(p, z - z)i2(Q) = f (5J_(P + C^p, 4>), 

Jo 

where we used the terminal and initial conditions p{T) = 0 and (()(0) = 0, respectively, 
which are well defined in view of Remark 2.6 On the other hand, setting (j) as test 
function in the weak version of (3.8) and integrating in time yields 

rT 


/ {d^_tP +C^P,(I)) = (U- Ud,(l))L^Q)- 
Jo 


The desired identity (3.11) follows easily from the derived expressions. □ 

We now prove necessary and sufficient optimality conditions for (1.1)-(1.3). 
Theorem 3.9 (first-order optimality conditions), z G Zad is the optimal control 
of problem (fllll-dOl) if and only if it solves ( |3.9| ), where p = p(z) solves ( |3.8[ ). 

Proof. We recall the control to state operator S : L^(0, T, IHI“®(D)) — >■ W defined 
by S(z) = u(z), where u(z) G W solves problem ( 1.2 ). Next, we write S(z) = So(z)-l-V^O: 
where So(z) denotes the solution to (1.2) with f = 0 and Uq = 0, while tpo solves 
with z = 0. Since Sq is linear, in our setting the variational inequality (3.10) reads 

(U - So(z - z))i 2 (Q) -I- fl{z, Z - z) 2 , 2 (q) > 0. 

Since So(z — z) = Sgz + il^o — {ifo + Sqz) = u(z) — Q, the previous expression becomes 

(0 - Ud, U(Z) - 0)i2(Q) -P P.{Z,Z- z)i2(Q) > 0. 


Using identity (3.11) of Lemma |3.8[ we arrive at 

(P, z - z)l 2 (q) -P /r(z, z - z)i 2 (Q) > 0, 


which is (3.9) and concludes the proof. □ 

3.3. Regularity of the optimal control. Since we shall be concerned with 
approximating the solution to the control problem (1.1)-(1.3|, it is essential to study 
its regularity. Here, on the basis of a bootstrap argument, we obtain such results. 

In what follows we will, without explicit mention, make the following regularity 
assumption concerning the domain D: 


kllH2(n) ^ II'Cw||l 2 (q), Vw G iL^(D) n i?o(D), 


(3.13) 


which is valid, for instance, if the domain Q is convex m- In addition, we will need 
the following assumption on a and b defining the set Zad^ 


a < 0 < b on dfl x (0, T). 


(3.14) 
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Theorem 3.10 (regularity of z). Let 7 = 1 . Assume that, for every e > 0, we 
have ^,Ud G and that Ug G //a,b € H^{Q) and (3.14| 

holds then z, the solution to the optimal control problem (flriji-ifOt, satisfies 

PI|L2(0,T;Hi(n)) + ||z|lRi(0,T;L2(n)) ^ ||f||L2(0,T;Hi-7n)) 

+ llud||L2(o_r.jji-£(Q)) + ||uo|iHi-'(a)) 

where the hidden constant does not depend on the problem data. Moreover, z S 

Proof. The proof is based on a bootstrap argument as in [31 Lemma 4.9], so we 
merely sketch it. By assumption, the right hand side of the state equation (1.2) satis¬ 


fies f-|-z G Lf{Q), while the initial condition satisfies Ug G IHI^(n). Standard regularity 
arguments yield that the solution verifies u G iL^(0, T; L^(n)) n L°°(0,T;IHI®(n)). In 
addition, by writing the problem as 

T®u = f -I- z - 5tD G 


we realize that u G L^(0, T; IHI^®(n)). The right hand side of (3.8) verifies u — G 
L'^{Q) so p has the same regularity that u possesses. From [35l §3.6.3] we have 


z = max ■( a, min ■( b,-p 


(3.15) 


This immediately yields z G H^{Q,T-,L'^(VL)). 

To obtain the claimed space regularity we recall that p G L^(0,r;]HI^®(fl)) and 
consider two cases: 

s G [ 5 , 1 ): Since p G L^(0,T; iLg (11)), formula (3.15) yields z G L^(0, T; iLd(ll)). 
Notice that assumption (3.14) is needed here to preserve the boundary values. 
lAI s G (0, 5 ): We now begin the bootstrapping argument. A nonlinear operator 
interpolation argument as in m Lemma 4.9] yields that z G L^(0,T;iL^'®(ll)). From 


this and condition ( 3.14| ) we conclude that z G L^(0, T; IHI^'®(11)). Define W 3 = £®u and 
notice that, since f-l-z G L^(0, T; IHI^®(D)), W 3 solves dtW 3 + = £®(f-|-z) G L‘^{Q) 
with initial condition u' 3 ( 0 ) = £®ug G L^(D). Therefore W 3 G L^(0,T; ^(ll)), which 
implies u = £~^W 3 G L^(0, T; IHI^®(11)). Define now 93 = £®p and notice that the same 
arguments yield that 93 G L^(0, T; IHI®(11)) and, therefore p G L^(0,T;]HI^®(11)). We 
consi der, again, two cases: 

s G [|, |): As in step Q] we have that z G L^{0,T-, Hq{LI)). 
s G (0, i): Nonlinear operator interpolation and (3.14), again, give us that z G 


2.1 


2.2 


L2(0,T;H3«(D)). Define now W 4 . = which solves dtW 4 +C^W 4 , = £3«/2(f + z) g 

Lfi{Q) with r(;4(0) = £^®/^Ug G Lf{Vl). This again yields that u,p G L^(0,r;]HI'^®(ll)). 
We con sider, one more time, two cases: 

s G [j, |): In this case z G L^(0, T; iLd(^))- 
s G (0, |): Define W 3 = and argue as before. 


2 . 2.1 

2 . 2.2 


Proceeding in this way we can conclude, after a finite number of steps, that for 
any s G (0, |) we have z G L^(0, T; 17^(11) )■ This concludes the proof. □ 


Remark 3.11 (regularity of u and p). Notice that while proving Theorem 3.10 
we have also shown that u, p G L^(D)) n L^(0, T; Hq{LI)). 

3.4. The extended control problem. To circumvent the nonlocality of the 


operator in problem (1.1)-(1.3) we realize it using the Caffarelli-Silvestre extension. 
In what follows we consider the equivalent problem: Find min{ J(trf 2 z) : G 
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V, z S Zad} subject to the extended state equation: Find G V such that tro ^ (0) = 
Uq in and, for a.e. t G (0,T), 

(trna7‘^,tra(/))+ = (f+ z,tra(/)) Wcj) G ,C). (3.16) 


To describe the optimality conditions we introduce the extended adjoint problem: 
Find gN such that trn = 0 in 17 and, for a.e. t G (0, T), 

(trn9j_(.^,trn^) +a(^,j>) = (tro- u^, tro ^)l 2 (o), G Hl(y°‘,C). (3.17) 

The optimality conditions in this setting now read as follows: the pair ('^(z),z) G 
V X Zad is optimal if and only if ^ (z) solves (|3.16 ) and 


(tro ^ + piz,z- z)i 2 (Q) >0 Vz G Zad, 


(3.18) 


where = J^(z) G V solves (3.17). 


4. A truncated optimal control problem. The state equation (3.16) is posed 
on the inhnite cylinder C = fl x (0,oo), therefore it cannot be directly approximated 
with finite element-like techniques. The first step towards discretization is to truncate 
C to a bounded cylinder Cy = fix (0, ff), which is possible because the energy decreases 
exponentially in [y ; see EH Proposition 4.1] for details. 

Proposition 4.1 (exponential decay). If, for a given 7 G (0,1] and s G (0,1), 
^ {z) gN solves (3.16), then for every jy > 1 we have 


ll^'^llL 2 (o,r;L 2 (y“,ax( 7 .oo))) e f + z). 


(4.1) 


where is defined in (2.13). 


Proposition |4.1| motivates a truncated control problem as follows. We first define 

Hl{y‘^,Cy) = {w G H^{y°‘,Cy) : w = 0 on dyCy U fl x {(T}} , 

Yy = {wG L\Q,T-Hl{y^,Cy)) : df imw G 7.2(0, T; H-'*(fl))}, 


and, for w,4> G Hf{y°‘,Cy), the bilinear form 


ay{w,(f)) = ^ f y°‘{A{x')'^w ■'V(j) + c{x')w(j)) dx'dy. (4.2) 

“s Jcy 

We define the truncated control problem as: Find min{ J(trn v,r) : v G Yy, r G Zad}, 
subject to the truncated state equation: Find v G Yy with tro t(0) = Uq in fl and 

(tro57u,tro(/>)-ko^(u,(/)) = (f-|-r,tro<()), ^(j) G ,Cy), a.e. t G {0,T). (4.3) 

As an instrument we define Ha '■ IHI®(fl) — Hj^{y°‘,Cy), the a-harmonic extension 
to Cy, i.e., if w G IHI®(fl), then w = HaW solves 


div{y°‘Yw) = 0 in Cy, w = 0 on dyCy U fl x {7"}, w = w on fl x {0}. (4.4) 


Remark 4.2 (initial datum). The initial datum Ug of (1.2) determines u(0) only 
on fl X {0} in a trace sense. We thus define u(0) = "HqUo. Remark 3.4 in provides 


the estimate ||Vu(0)||L2(yo < Huolln^n)- 
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Let us now provide, for 7 = 1, an energy estimate that will be useful to derive an 


L^{Q) error estimate for the fully-discrete scheme of ^5.3 


Theorem 4.3 (energy estimate: 7 = 1 ). Let s S (0,1) and 7=1 and denote by 
V S Yy the solution to ( |4.3| ). Iff and r belong to L‘^{Q) and Uq S then 

II trn 9t'y||L2(Q) + ||Vu||i=o(o,T;L2(y-,c,r)) - 11^ + A\l^{Q) + IIuoIIh^o), (4.5) 

where the hidden eonstant does not depend on v nor the problem data. 

Proof. Set (f) = dtv in (4.3), integrate over time and use the estimate of Re¬ 
mark]^ |lV?;( 0 )||L 2 (y._c ,4 ^ IIuoIIh^o)- □ 

As m ( |3.2| we introduce the truneated adjoint problem: Find p G Yy such that 
trop(T) = 0 and, for a.e. t G (0,T), 

{tin Off _^p, tin (j) + ay{p,(j)) = {trnv - Ud,tTQ(j)), 'i(j) G Hy{y°‘,Cy)- (4-6) 

The same arguments provided in Theorem |3.9| allow us to conclude that the pair 


(4.7) 


(w(r), r) G Yy X Zad is optimal if and only if u(r) solves (4.3) and r satisfies 
(trnp + /rr, r-f)i 2 (Q) > 0 Vr € Zad, 


where p = p{r) G Yy solves (4.6). 

The next result shows how (h(r), r) approximates ('^(z),z). 
Lemma 4.4 (exponential convergence). For every ‘J >1 we have 


I*"- z||l2(Q) < e“^“^(A^(uo,f) + ||r||L2(0,T;M--(O)) + ||Utz||L2(0,T;M—(O))), 


and 


W'^ — ?;| 


L2(o.T;Ai(y“,C)) ~ ® ^ ^(^7(^0, 0 

+ l|f||L2(o,T;H-7n)) + l|Ud||L2(o,T;H-7n)))) 


(4.8) 


(4.9) 


where A.y is defined in (2.13). 

Proof. We proceed in four steps: 

[T] Set z = r G Zad and r = z G Zad in the variational inequalities (3.18) and (4.7), 
respectively, and add the obtained inequalities to arrive at 

Ml|r-z||i2(Q) < (tro(#-p),r-z)i2(Q) 

= (tro(.# - .^(r)), r - z)i 2 (Q) + (trni^ir) - p),f - z)i 2 (Q). 

lAI Consider (tro(# — .^(r)), r — z)y 2 ^Qy Define fj := tP — {?) G Y and observe 

that tro V'(^) = 0 and, for all (fp G Hy{y°‘,C), we have 

J + a{tp,4'p)) dt = J -'P{T)),tTn(j)p)dt. 

Analogously, define (p := "P — ^{r) gY, which satisfies tro v?(0) = 0 and 

r r - - 

/ ((trn 97 ^ 3 , trn </>«) + 0 ( 75 , </>«)) df = / (z - r, trn (('„)i 2 (n) dt, V(/)„ G 77i(y“,C). 

Set 4>p = 4>u = and apply Lemma 2.7 Since tin 'ip{T) = tr^ (/?(0) = 0 we get 

(tiniP - P(J)),7 -z)L2f^Q) =-||tra(‘^-'i'(r))|l|2(Q) <0. 
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lAI The previous step shows that proving (4.8 1 reduces to obtaining a bound for 
II ~p)||^ 2 (Q) and [571 Lemma 4.3] yields such a bound. 

lAI In order to prove ( |4.9| ) we write 

'^(z) — i;(r) = (^(z) — ^ (r)) + (fll (r) — i;(r)). 


The first term satisfies (2.11) with right hand side z — r so that by (4.8) this term is 


bounded. For the second term we again apply [571 Le mma 4.3]. □ 

Remark 4.5 (regularity of r vs. z). In Theorem 3.10 we studied the regularity 
of z. The techniques of EH Remark 4.4] allow us to transfer these results to r, 
the solution of the truncated optimal control problem. In a similar fashion, we can 
establish the regularity results of Remark 3.11[ for trn v and tr^p. For brevity we skip 
the details. 

5. Approximation of the state equation. We recall the numerical approxi¬ 


mation of the state equation (2.11) developed in [57]. The scheme employs first degree 
tensor product finite elements in space and finite differences in time. The latter is the 
backward Euler scheme for 7 = 1 whereas, for 7 G (0,1), it is the scheme of [50l [5T| . 
which was studied under appropriate time-regularity conditions on the solution ^ 
in m We also derive a novel L^{Q) a priori error estimate for the fully discrete 
approximation of the state equation (2.11) with 7 = 1 and s G ( 0 ,1). 


5.1. Time discretization. Let /C G N denote the number of time steps. Define 
the uniform time step r = T/K, > 0, and set tk = kr for Q < k < K,. We denote 
the time partition by T := {tk}k=o- H £ C{[Q,T],X), we denote (j) = (j){tk) and 
(jf = {<(>^}^_g. On such sequences we define the norms 


w 


.(x) = max Wcj) 
^ ’ o<fe</c 


YZ 


■u' 


1 ^- 


Over sequences (jX C X we define the discrete time derivative by 

fc = 0,...,^-l. (5.1) 

As in m §3.2] we also define, for 7 G (0,1), the discrete fractional derivative as 


F(2 - 

1=0 


fc +1 k — 1 

i^k+1-3 _ z _^ Qj - Qj+i ^fc-i _ 

1=0 




(5.2) 


where Oj = (j -I- 1 )^“''' _ j^~'y and provided the sum for A: = 0 is defined to be zero. 

We remark that, any sequence (f>'^ G X can be equivalently understood as a 
piecewise constant, in time, function (j) G L°° {0,T; X): 

= (j)’", VtG (4-1,4], fc = l,...,/C. (5.3) 


This identification will be very useful and, in what follows, we will use it repeatedly 
and without explicit mention. 

5.2. Space discretization. The space discretization is based on truncation and 
the finite element method. The truncation is as in m Lemma 4.3], which shows that 
truncating C to Cy induces an exponentially small error. Since we are now dealing 
with the bounded domain Cy^ we can discretize using finite elements. 
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The finite element discretization follows [29l §4]. Let = {K} be a conforming 
triangulation of fl into cells K (simplices or n-rectangles). We denote by Tq the 
collection of all conforming refinements of an original mesh and assume Tq is 
shape regular [S]. If S Tq we define hj^- We define .5^ to be a 

partition of Cy into cells of the form T = K x I, where K G S/q, , and I is an interval 
that comes from the partition {ym}m=o defined by 


Um — 




m = 0 . M, 


(5.4) 


where ( = (((a) > 3/(1 — a) > 1. The set of all such triangulations is denoted by T. 
Note that the following weak regularity condition is valid: there is a constant a such 
that, for all G T, if Ti = iLi x /i, T 2 = K 2 x I 2 G 3'y have nonempty intersection, 
then < a, where hi = |/|; see [551 [H]. For ,5^ G T, we denote by 0\i{^y) the 

set of its nodes and ^{S^y) the set of its interior and Neumann nodes, respectively. 
We also denote by = #5^(^) the number of degrees of freedom of ^y. We assume 
that #,% « M” so that N « M”+b 


The main motivation to consider elements as in (5.4) is to compensate the rather 
singular behavior of solution to problem (2.11) as y « 0+; see [57] for details. 


For h7y G T and F^ = dLCy U 11 x {(T} we define the finite element space 

N{^y) = {W G C{Cy) : W\t G Vi{K)®f‘i[I) \/T = K x I G ^Ty, lF|r^ = O} , 

If iF is a simplex, then 'Pi{K) = Pi(iL), whereas if iF is a cube, then ViiK) = Qi{K). 
We also define U(c%) = trnY{£^y), i.e., a Vi finite element space over the mesh 


5.3. A fully discrete scheme. The fully discrete scheme to solve ( 1.6 ) combines 
the space discretization of §5.2| with the time discretization of §5.1[ To define it, we 
first consider the weighted elliptic projector studied in §4.3]: 

w G ,Cy) : ay(Gs^w,W) = ay{w,W), VW G¥(^y). (5.5) 

The fully-discrete scheme computes C Y{3^y), an approximation of the solu¬ 


tion to (4.3), with r = 0, at each time step. We initialize the scheme by setting 

O TIqUo, 


(5.6) 


where ° hia and T-La is the a-harmonic extension operator defined in 

notice that tro = tro G^yU(O), where r!(0) solves (4.4) with w = Uq. 

For fc = 0,... ,/C — 1, G V(/5y) solves 

{5'^ivnV^+\tvnW)L.(^n)+ay{V^+\W) = {‘i^+\tTaW), VW G ¥(.%), (5.7) 

where 5"* is defined by ( |5.2| ) for 7 G (0, 1) and by ( |5.1| ) for 7 = 1 and f*^+^ = 
T~^ f dt. An approximate solution of problem (1.2) is C U(C?o) with 




(5.8) 


To present error estimates for scheme (5.6)-(5.7), for 7 G (0,1), we introduce 


A = A(v,g) = |lv||H»(a) -I- llg|l//2(o_T;H-»(a)), 

B = B(y,g) = ||v||Hi+3«(n) -I- |lg|t=o||Hi+'>(n) + llsll 
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where z/ > 0 is arbitrary. Theorem 5.3 in 122 ! provides the following error estimates 
for the scheme (|5.6|)-(5.7 1 with 7 € (0,1) and s € (0,1). 


Theorem 5.1 (error estimates: 5,7 G (0,1)). Let r = 0 and 3,7 G (0,1). If v 


solves (4.3|, VJ solves (5.6)-(5.7), ^(uo,f)S(uo,f) < 00 and SLy verifies (|5.4[), then 


|2 ^ I /\H2s 1 —— 


tra(^;" - (T)] 2 < t^A{uo, f) + \ log 7V|^«iV—^(uo, /), (5.9) 


and 


11^" - ^ ^ f) + I log iVriV;l+Te(uo, /), 


(5.10) 


where 0 G (0, ^), and the hidden constant does not depend on v, Vy nor the problem 
data, but blows up as 9 f 

5.4. L^{Q)-error estimate: s G (0,1) and 7 = 1 . We now derive a novel 
L^(Q)-error estimate for (5.6)-( [5)7| ) with 7 = 1 , which is inspired by classical tech¬ 
niques developed, for instance, in [illlS]. To obtain it, we set r = 0 and consider, as 
a technical instrument, a semi-discrete approximation to (4.3): Set = T-La^o and, 
for fc = 0,...,/C — 1, compute 1/^+^ G Hf{y°‘,Cy), that solves 

(54rnl^''+\trn(())L2(o) + ay{V'"+'^,(j)) = G Hl{y^,ey), (5.11) 

where is defined by We present the following stability result. 

Lemma 5.2 (stability). Let C Hj^{y°‘,Cy) solve (5.11). //f G Lfi{Q) and 
Uq G IHI®(n), then we have 

II tro (5^W||^ 2 (i 2 (Q)) -I- IIVy^||^oo(i 2 (j^a^Cy)) < ||f||^ 2 (i 2 (n)) -I- ||uo||Ho(n)) (5-12) 

where the hidden constant does not depend on nor the pro blem data. 

Proof. Set (f = 14^+^ — and use the estimate of Remark |4.2| □ 

Define the piecewise linear function V G {[0,T]] Hf{y°‘,Cy)) by V{0) = 


v{t) = v>^ + it- tk)S^V’^+\ t G ( 4 , 4 + 1 ], 


(5.13) 


for fc = 0,...,/C — 1. Using this notation, we rewrite equation (5.11) as 


(traatU(t),trn((i)i 2 (o) + ay{V^{t),(j)) = (r(t),tra 0), yfi G Hl{y°‘,Cy), (5.14) 


for a.e. t G (0,T). We are now in position to derive an error estimate for (5.11). 

Theorem 5.3 (semi-discrete error estimate: 7 = 1 ). Let v and solve (4.3) 
and (5.11), respectively. // f G L°°(0, T; L^(D)) and Uq G IHI^(r2), then 

II trn(ii - W)||L2(o,T;L2(n)) < t (||f||L“(o,T;L2(n)) + lluollH^n)) 5 (5.15) 

where the hidden constant does not depend on v, nor the problem data. 

Proof. Define e = v — V and e = v — V^. Set 7 = 1 and r = 0 in (4.3) and then 
subtract from it (5.14). Integrating with respect to time the result we obtain 


(tro e(t), tro 4>)L^(n) + ay 


e(?)dC,(() = 


-f (trn(e(<) -e(f)),tra fi) i^q.), 


a.e. tG(0,T), yfi G Al{y°,Cy). 
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Set, for a.e. t G (0,T), cj) = e{t) G Integrating over time once more yields 


|trae(t)||i 2 (n) dt < 




(trn(e(t) - e(t)), trn e(t))i 2 (Q) dt 


(5.16) 


where we used that 


J ay(^Jey)d^,e{t)^dt=^a.y^J e{t)dt, e(t) > 0 . 

Notice that, since = t~^ 

r im - de=E r - f") 

*^0 k—l 


Consequently, if t; < t < tj+i, we have 


[\m - f"(C) d? = fm) - f"(e)) d? < r||fiu„(o.*). 

Jo Jti 


In conclusion, the first term on the right hand side of (5.16) can be bound by 

rT / rt 




li“(o,T,L 2 (n)) + ^11 trn e||i 2 (Q). 


Since, on (tfe,tfe+i], we have that |e(t) — e(t)| < estimate (5.12) yields 

^ II trn(e(t) - e{t))\\l-,^^^ dt < {\\^^\\%(LHn)) + lluollH-(n)) ^ 

and, therefore, the second term on the right hand side of (5.16) can be bounded by 

rT 


'0 


(trn(e(t) - e(t)), trn e{t))L 2 ^Q) dt 


< ^11 trn e"^ 


L2(o,T;L2(n)) 

+ (l|r'll^2(L2(n)) + lluollHocn)) ■ 


Collecting all the derived bounds we arrive at the desired error estimate (5.15). □ 
With this estimate at hand we can control the difference between the fully and 
the semi-discrete problems. 

Theorem 5.4 (auxiliary error estimate: 7 = 1 ). Let 7=1 and assume that 
uo G Hi+®(n) and f G L^{0,T;B.^-%n)). If and solve ( [KTI| ) and ([^, 
respectively, then 


trn(14'^ — VJy)||£2(j;,2(n)) < |logA^p®iV "+i (||C||^2(o_T;Hi-=(n)) + lluollni+^n)) j 


where the hidden constant does not depend on v, nor the problem data. 
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Proof. We start by defining the error 


= (V^ - - VL) = r + . 


(5.17) 


where is defined in (5.5). We estimate 9'^ by invoking the approximation prop¬ 


erties EH Proposition 4.7] of and the regularity results of EH Theorem 2.7]: 

II tra 6i'^||£2(i2(n)) < | logfVp®fV“^^+i (||P||^2(]ai-»(o)) -I- ||uo||Hi+'>(n)), 

The estimate of follows along the same lines of El Lemma 5.6]. For brevity, 
we skip the details. □ 

We collect the estimates of Theorems 5.3 and 5.4 to derive a {Q)-eiror estimate. 


Theorem 5.5 (error estimate for v: 7 = 1). Assume that 7 = 1 and let v and 
VZf solve (4.3) and (5.7), respectively. If ^ G L'^(0,T;L^(r2)) n L^( 0 ,T;IHI^“'*(n)) 


and ug S IHI^“''®(n), then 

||tro(r;-V'^^)||i2(Q) <t 


L°°( 0 .r;L 2 (o)) -I- ||uo||Hs(n)) 


+ |logfVp^A^ "+1 (||f||L2(0 7'.Hi-.(n)) + ||Uo||hi+<.(o)) , 


where the hidden constant does not depend on v, nor the problem data. 

Corollary 5.6 (error estimate for u: 7 = 1 ). Assume that 7 = 1 and let 
u solve ( |1.2[ ) with z = 0 and be defined by (5.8). If f G L°°(0, T; L^(f2)) n 
L^(0, T; EI^“'’(ri)) and Ug e IHI^+®(n), then 


1 ^ “ t^5hllL^(0.T;L2(n)) 


< 


U“(0,T;L2(O)) -I- ||Ug||Hs(o)) 

lOgfVP^lV " + 1 (||f||L 2 (g_y.]gl-.(f 2 )) -I- ||Ug||]gil + s(Q)) , 


where the hidden constant does not depend on v, nor the problem data. 

Proof. The result is a consequence of Theorem |5.5| and EH Lemma 4.3], in 
conjunction with the appropriate choice of fh explored in |29L Remark 5.5]. □ 

6. Approximation of the fractional control problem. We propose an im¬ 
plicit fully-discrete scheme to approximate the solution of the fractional control prob¬ 


lem (1.1)-(1.3): piecewise constant functions for the control and, for the state, first 
degree tensor product finite elements in space, as described in §5.2[ and the finite 
difference discretization in time detailed in 95.11 


As stated in Theorem 5.1 in order to have the error estimates (5.9) and (5.10) for 


the approximation of the state equation ( 1 . 2 ), we have to require that A(uo, f-|-r) < 00 . 
This strong in time regularity assumption is not satisfied by the optimal control 
r, meaning that we are not able to apply the results of Theorem |5.1[ This is in 


sharp contrast with the case 7 = 1 which, according to Theorem 5.5 only requires 
f,f G L°°(0,T; n L^(0, T; IHI^“®(n)) which, by imposing (3.14) and invoking 

Remark |4.5| and Theorem |3.10[ is satisfied by the optimal control f. Due to this 
regularity restriction we can obtain an error analysis for 7 = 1 only. We remark that 
L^(Q)-eiwr estimates, for s,'y G (0,1) are not available in the literature, especially 
under the correct regularity assumptions. In §6.2| we will present error estimates for 


7=1 and s G (0,1), and in (6.3 we will show the convergence, without rates, for the 


remaining range of parameters. 

Finally, to simplify the exposition, in what follows we assume that a and b are 


constants that satisfy (3.14). 
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6.1. An implicit fully discrete-scheme. To discretize the control we intro¬ 
duce the finite element space of piecewise constant functions over 

= {Zg L°^{n) : Z\k G Po(A), VK g , 

and the space of piecewise constant functions in time and space 

Z(r, = {z^ c L°°(Q) : G Z(5h)} ■ (6.1) 

We define the space of discrete admissible controls as follows: 

Zad(r,.%) = Zadnz(r,.%). (6.2) 

It will be useful to introduce the L^(Q)-orthogonal projection onto Z{T,^q). 
The operator IlJ-^ : L‘^{Q) —t Z(T, ^n) is defined by 


r G L^iQ) : (r - =0 VZ G Z(r, 

and, for all r G il^(0, T; L^(r2)) n L^(0, T; iJ^(n)), satisfies: 

W^-^So^WlHQ) ^ hsnW'^x'rWL^Q) +T\\dtr\\L2(^Qy 

Notice also that, since a and b are constant, IlJ-^Zad C Zad(T, .%). 
We define the discrete functional 


(6.3) 


(6.4) 


1 


JZ— -|| trn — u5||f2(i2 


m) 


+ ^\\Z^ "2 




where the .^^-norm is defined in (5.1 The identification between a sequence (j)'^ and 
the piecewise constant function (5.3) will be used repeatedly below. For instance, if 
uj = Ud we would have that J'^^{w,r) = J{w,r) whenever w'^ = w and = r, that 
is, the arguments are piecewise constant over T. This was already implicitly used in 
( 6 . 2 ), when we defined Zad(T, ^n). 

The numerical scheme reads: Find min subject to the discrete 

state equation: initialize as in (5.6) and for fc = 0,...,/C — 1, let G V(i3^) solve 


(^"^ tro V^+\trn M^)L=(n) + a^iV^+\W) = (f'^+i + tm W) 


(6.5) 


for all W G V(5y) and the control constraints Z'^^^ C Zad(r,.%). If(y^^,Z5.J 
denote the solution to this problem, setting 


= trn w 


.5V 


( 6 . 6 ) 


we obtain a fully-discrete approximation {U^^, Z'^^) G x Zad(T,.%) to the 

fractional control problem (flT|)-([L3t. 

Remark 6.1 (locality). The main advantage of the scheme (6.5) approximating 
the fractional control problem ([Tt|)-([L31) via ( |6.6[ ) is its local nature. 

6.2. A priori error analysis: 7=1 and s G (0,1). Let us consider s G (0,1) 
and 7 = 1 in ( 3.16)-(3.17) and provide an a priori error analysis for the fully-discrete 
scheme proposed in §6.1| To do so, we provide first order necessary and sufficient 
optimality conditions of the fully-discrete problem. We define the discrete adjoint 
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problem: Find C V(^) such that trn P^^ = 0, and for k = JC — 1,... ,0, 
G solves 

{SHin P^^,tTQW)L2(n) + ay{P^^,W) = - u^+\tro , (6.7) 


for all IF G 
condition reads: (F^^, 


Here denotes = —t ^ (<)) 

^ 5 h) optimal if and only if FJ^ 


. The optimality 
solves (|5.6| and (|6.5|) and 


{ti'Q P^^ + Z — Z'^^)[^2(^Q'f > 0 VZ G Zad(T, =%), ( 6 . 8 ) 

where P^^ solves ( |6.7[ ). Notice that ( | 6 . 8 | ) can be equivalently written as 

(trnP^^ + A^Z^,,Z-Z|.Ji.(n) VZ G Z(.%), a < Z < b, Vfc = l,...,/C. 

To see this, it suffices to set Z'^ = Zx(tk-i,tk], with Z G Z(,%) and a < Z < b. This 
greatly simplifies the implementation. 

Let us now introduce two auxiliary problems. The first one reads: Find C 

V{^y) such that trn = 0 and, for fc = Ail — 1,..., 0, G V(=i5y) solves 


((5^ tro bF)i 2 (f 2 ) + ay{Q%^,W) = (tro - u^+\ tro W) , (6.9) 

for all IF G V(^), and where v = u(r) solves ( |4.3[ ). The second one is: Find 
C Y{3^y) such that tr^ = 0, and for fc = /C — 1,..., 0, G V(5^) solves 


((54raP^^,trnTF)i2(Q) +a^(P|.^,fF) = (|tra F^+^(r) - u^+\trnfF^, (6.10) 

for all IF G Y{^y). These auxiliary problems will allow us to derive error estimates 
for the fully-discrete scheme proposed in § 6 . 1 [ 

Lemma 6.2 (error estimate for the control: 7=1 and s G (0,1)). Let r be the 
solution to the truncated optimal control problem of and let ZJ.^ be the solution 
to the fully-discrete optimal control problem of ) 6.1 Assume that Uq G IHI^+'’(f2) and, 
for every e > 0, G L^{G,T■,W^-\n))C^ H^{Q,T-L'^{n)) ond f G . 

Then 

||r-Z5.JU2(Q)<r+|log7V|2^1V-^, 

where the hidden eonstant is independent of the discretization parameters but depends 
on the problem data. 

Proof. We proceed in several steps. 

|T] Setting r = ZJ.^ and Z = H^-^r in ( |4.7[ ) and ( | 6 . 8 [ ), respectively, and adding the 
derived inequalities we arrive at 

Mlir- ^5hlli2(Q) ^ (tra(p- + (FnP^y + r)i 2 (Q), 

where HJ^ is defined in 

Using the solution to (6.9| wewritep —P^^ = (p—) + (Q5y ~The first 
term is estimated by using the results of Theorem |5.5| as follows: 


II trn(p - Q 5 y)||L 2 (Q) < r (|| trn u||loo(o_t;L2(0)) + IIUd||Loo(o,T;L2(n)) + ||uo||H 7 n)) 
+ I logiVp®A^ "+i (II trn ll||L2(o.T;Hi-'’(n)) + l|ud||L2(o,T;Hi-=(n)) + l|uo||Hi+=(n)) ■ 
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Where we used that tr^ u € H L^(0,T;IH[^“®(il)), which follows from 

Remark 13.111 and Remark 14.51 

[ 3 ] To estimate the difference tro((5^^ — we write + 

— P'^ ), where solves ( |6.10 ). Employing the stability estimate established 
in |271 Lemma 5.1] and Theorem] 5.5 [ we arrive at 


II tra( 0 ^^ - R\)\\l^{q) < II trn(i; - r^^^(r))||L 2 (Q) 

^ T (||r|U“(0.T;L2(o)) + ||uo||H'>(n)) 

+ I logiVp®lV "+i (||r||i 2 (Q_y.jji-s(Qp + ||uo||Hi+«(n)) ■ 

To handle the term P^^ — P^^ we invoke the discrete counterpart of Step 2 in 
Lemma |4.4[ that is an argument based on summation by parts, to arrive at 


(tra(P5r^ — P^y), — P)l'^(q) < 0. 


[ 4 ^ Using the solutions to (6.91 and (6.10) we write 


(traPJ-^ r)i 2 (Q) = (trnp +/rr, - r)i 2 (Q) 

+ (trn(P^^ ± Q\ - P), nj-^r - r)i 2 (Q) + - r, - r)i 2 (Q) = I + II + III. 

Using the properties of the projector II]^^ and the smoothness of p and r we have 

1= (trap + pr-n]^^(trap + pr),nP^r-r) < (r|| tr^p + pr||^i(o,T;L2(n)) 


The term II can be handled by repeating the arguments of Steps 2 and 3, while III is 
controlled by a trivial aplicaiton of the Cauchy Schwarz inequality. 

[~^ The assertion follows from collecting all the estimates we obtained in previous 


steps and recalling that ^ N~ . □ 

On the basis of of Lemma [6.2] we derive the following important result. 
Theorem 6.3 (control error estimates: s G (0,1) and 7 = 1). Let z be the 
solution to the space-time fractional optimal control problem (!.!)-(1.3) and let 
be the solution to the fully-discrete optimal control problem o/lj In the framework 
of Lemma \ 6. we have the following error estimate 


llz-^5-JU2(Q)<r+|logfVp«iV-^ 


where the hidden constant is independent of the discretization parameters hut depends 
on the problem data. 

Proof. The result follows from Lemmas |4.4| and |6.2] in conjunction with an 
appropriate selection of the parameter ZL, as detailed in Remark 5.5] and [21 
Corollary 5.17]. □ 

We conclude with an error estimate for the state in the L^(0,T;]HI®(O))-norm. 

Theorem 6.4 (state error estimates: s G (0,1) and 7 = 1 ). Let u be the optimal 
state of the space-time fractional optimal control problem (1.1)-(1.3| and let be 
defined as in ( |6.6| ). In the framework of Lemma 6.2 we have the following error 
estimate 


l|u — llL2(0,T;H'’(n)) + I log A^P®iV"+i 
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where the hidden constant is independent of the discretization parameters hut depends 
on the problem data. 

Proof. We first write 


|u - Usn\\L^{0,T-U‘’{Q)) < || - 'y)||L2(o,T;H=(n)) + || tro V - C/5f;||L2(o_r;Hs(f2)), 


and note that the first term is controlled in (4.9). The second term is handled by 
noticing that 


II trn w — {75-^11^2(0,T;H'’(a)) — II tra(r; — V|5K)lli=(o,T;H»(n)) 

< II trn(7 - i54-(''))IU=(o,T;H=(n)) + II “ K^y)l|L2(o.T;H=(n))- 

Since f+ f G 77^(0, T; L^(f2)), the first term on the right hand side of this inequality is 
estimated using the error estimates for the discrete scheme presented in EH Theorem 
5.4]. The second one can be handled by invoking the stability of the discrete scheme 
and the error estimates of Theorem |6.3| Collecting these bounds we obtain the result. 
□ 


6.3. Convergence. Let us now consider the case when either 7 , s G (0,1), or 
the problem data is not smooth enough to yield the error estimates of § 6 . 2 | and elu¬ 
cidate the general convergence properties of the fully discrete scheme. Notice that 
we are not only approximating the state equation via discretization, but we are also 
approximating the cost, so convergence of discrete optimal controls to the continuous 
one is not immediate. To begin, as in Definition |3.1[ we introduce the discrete control 
to state operator 


: Z(r, 5 h) 9 ^ VJ, C Ym, 


( 6 . 11 ) 


where VJy solves the discrete state equation (6.5). Notice that the stability estimates 
implicit in Theorems 5.1 and 5.3 \m Lemma 5.1]) yield that, for all 7 , s G (0,1) 


the family {S]5y} is uniformly bounded for f, ^ Moreover, the error 

estimates imply the pointwise convergence of these operators so that, by the uniform 
boundedness principle, they converge uniformly to S. This will be crucial in showing 
convergence. 

With the discrete control to state operators at hand, like in (3.2), we define the 
reduced cost functional by 


fUz. 


% > 


Jl isl z 




1 , 


= Alltrn - u 


Jdll^2(L2(n)) 


All^5hllt2(L2(n)). 


( 6 . 12 ) 


The convergence of the fully discrete scheme is the content of the next result. 
Theorem 6.5 (convergence). The family 'Is uniformly bounded 

and it contains a subsequence that converges Lf{Q)-weak to r, the solution to the 
truncated optimal control problem. Moreover, 1 / 7=1 the convergence is strong. 

Proof. Boundedness follows immediately from the fact that minimizes 
If zq G Zad, then 


F'^y.{Z'fra) ^ ^ lko|li2(Q) + ||ud||i2(Q), 

where we used the uniform boundedness of and Sp^. This implies the existence 
of a (not relabeled) weakly convergent subsequence. 







A fractional space-time optimal control problem 


23 


To show convergence of this subsequence to r we appeal to the theory of T- 
convergence, for which we need to verify several assumptions: 

|T] Lower hound inequality. Assume that ^ 2 in L'^{Q). For w G L‘^{Q), we have 

- S2,u;)i2(Q) = Sz,u;)i2(Q) + - z ), w ) l ^ q ) = I + II. 

The pointwise convergence of to S shows that I —>■ 0, while their uniform conver¬ 
gence that II -G 0. In conclusion Sz. Lower semicontinuity of the norms 

and uj —>• \id in L‘^{Q) imply 

f{z) < liminfFj^(Z 5 .^), 
which is what we needed to show. 

[~^ Existence of a recovery sequence: Let z G Zad, then nj^z G Zad(T,.%) con¬ 
verges strongly to z in L‘^{Q). Consequently, Sj^nj^z —>■ Sz in L‘^{Q) as well. The 
continuity of then implies 

/(z) > limsupFj^(n 5 -^z). 


[~3~| Equicoerciveness: Since 


we have, by [9l Proposition 7.7], that the family {F^} is equicoercive. 

[ 4 ] Steps 1 and 2 show the T-convergence of the discrete reduced costs F]^^ to the 
reduced cost /. This implies, using [HI Corollary 7.20], that minimizers of F^, if they 
converge, must do so to a minimizer of /. Step 3 and the uniqueness of the minimizer 
of the reduced cost / are the conditions for the fundamental lemma of T-convergence 
[HI Corollary 7.24]. In conclusion, converges weakly to f, the minimum of the 

truncated cost functional. 

We conclude with the strong convergence for the case of 7 = 1, which follows from 
the a priori estimates. Namely, a basic energy estimate for scheme (6.51, together with 
a slight modification of Lemma |5.2| imply that 

T - ^TS57^]%)llL“(T,T;L2(n)) + II Irn 11^2(0 < 1, 

where St indicates the backward shift (in time) operator. These are, according to [TTl 
Theorem 1], sufficient conditions for relative compactness of the family 

{ } ■% GTr! ,T>0. 


By passing to a subsequence we get strong convergence. Denote now by F^^ the 
solution to the discrete adjoint equations ( |6.7| ) with right hand side — uj. 


Stability of (6.7) implies that the sequence {F^^} converges strongly in Lf{Q). 

Set r = and Z = in ( |4.7| ) and ( | 6 . 8 | ), respectively. Adding the derived 
inequalities we arrive at 


~ ^5hlli,2(Q) ^ ~ + (trn FJ.^ +,IlJ^r - r)j^2(Q) 


where p and F^^ solve ( |4.6| ) and 
we immediately have that (tro F^^ + — 1 )^, 2 (q) —)• 0. We handle the 

remainding term using the strong convergence of F^^. 

This concludes the proof. □ 


(6.7), respectively. Using the definition of IlJ^ 
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Fig. 7.1. Computational rates of convergence for the control and state on anisotropic meshes 
for n = 2 and s = 0.2, 0.4, 0.6 and s = 0.8. For a fixed number of time steps, K = 1400, the left 
panel shows the decrease of the L^{0,T-, L^{Q,))-control error with respect to N and the right one 
that of the L^{0,T; L^{fl))-state error. In the control case we recover the rate For the state 

we observe a rate of . 


7. Numerical experiments. Let us illustrate the performance of the fully dis¬ 
crete scheme proposed in §6.1|for 7 = 1 and the error estimates derived in §6.2| 


7.1. Implementation. The implementation has been carried out in MATLAB®. 


The stiffness and mass matrices of the discrete system (6.5) are assembled exactly, and 


the respective forcing boundary term are computed by a quadrature formula which 
is exact for polynomials of degree 4. The resulting linear system is solved by using 
the built-in direct solver of MATLAB®. To solve the minimization problem, we use 
the projected BFGS method with Armijo line search; see m- The optimization al¬ 
gorithm is terminated when the ^^-norm of the projected gradient is less or equal to 
10 -®. 


To illustrate the error estimates of j 6.2 we need an exact solution to the fractional 
control problem ( |1.1[ )-( [L3| ). Let n = 2, ft = 1, ft = (0,1)^, and £ = —A. In this 
setting, the eigenpairs of £ are: 

Afe,; = 7r^(/c^-I- ipk,i{x[,X2) = sm{kTrx[)sm{l7rX2) fc, Z G N. 

Set u = e* sin(27ra:j) sin(27ra:2), which yields f = (1 -I- A| 3 ) 6 * sin(27ra;j) sin(27rx2) — 


Set also p = —/i(T — t)e‘sin(27rx'i) sin(27rx2). Dehnition 3.4 then yields = [l — 
ft{—l + (1 — 2 )(T" ~ t) }] e* s in(27rx2) sin(27rx2). Finally, we set a = 0 and b = 0.5. 

The projection formula (3.15) gives the value of z. This defines, for any s G (0,1), 


the data and solution to the optimal control problem ([lT|)-((L31). 

7.2. Convergence rates in space. Let K. = 1400. The asymptotic relations 




\u-U% 


.%lli2(0.T;L2(n)) 


N- 


are shown in Figure [7T] The left panel illustrates the quasi-optimal rate of conver¬ 
gence for the optimal control with respect to the number of degrees of freedom N 
for all choices of the parameter s considered. As noted in 13 EH, in order to recover 
optimality, the state and adjoint equations must be discretized with the anisotropic 


rehnement, in the extended dimension, dictated by (5.4). 


From Figure [7Tl we can also observe that the approximate optimal state converges 
with a rate N~ 3 . This rate is not discussed in this paper and will be part of a future 
work. The theoretical rate of convergence for the approximation of the optimal state 
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l|2-2i|l2(0,T;i2(n)) 



Fig. 7.2. Computational rates of convergence for the control on anisotropic meshes for n = 2 
and s = 0.2, 0.4, 0.6 and s = 0.8. For a fixed number degrees of freedom in space, N = 927828, 
the figure shows the decrease of the L^{0, T; (Q))-control error with respect to K. In this case we 
recover the rate K~^. 


is dictated by the results Corollary 


5.6 


N 3 ^ which, in fact, is a consequence of 


the error estimate derived in [571 Proposition 4.7]. 


7.3. Convergence rates in time. Let N = 927828. The asymptotic relation 

l|z - llL2(0,T;L2(n)) « 

is shown in Figure [7^ and illustrates the optimal decay rate in the control with respect 
to /C, for all choices of the parameter s considered. 
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